A Fully Integrated Transport Approach to Heavy Ion Reactions with an Intermediate 

Hydrodynamic Stage 
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We present a coupled Boltzmann and hydrodynamics approach to relativistic heavy ion reactions. 
This hybrid approach is based on the Ultra-relativistic Quantum Molecular Dynamics (UrQMD) 
transport approach with an intermediate hydrodynamical evolution for the hot and dense stage of the 
collision. Event-by-event fluctuations are directly taken into account via the non-equilibrium initial 
conditions generated by the initial collisions and string fragmentations in the microscopic UrQMD 
model. After a (3+l)-dimensional ideal hydrodynamic evolution, the hydrodynamical fields are 
mapped to hadrons via the Cooper-Frye equation and the subsequent hadronic cascade calculation 
within UrQMD proceeds to incorporate the important final state effects for a realistic freeze-out. 
This implementation allows to compare pure microscopic transport calculations with hydrodynamic 
calculations using exactly the same initial conditions and freeze-out procedure. The effects of the 
change in the underlying dynamics - ideal fluid dynamics vs. non-equilibrium transport theory - will 
be explored. The freeze-out and initial state parameter dependences are investigated for different 
observables. Furthermore, the time evolution of the baryon density and particle yields are discussed. 
We find that the final pion and proton multiplicities are lower in the hybrid model calculation due 
to the isentropic hydrodynamic expansion while the yields for strange particles are enhanced due 
to the local equilibrium in the hydrodynamic evolution. The results of the different calculations 
for the mean transverse mass excitation function, rapidity and transverse mass spectra for different 
particle species at three different beam energies are discussed in the context of the available data. 

PACS numbers: 25.75.-q, 24.10.Lx, 24.10.Nz, 25.75.Ag 
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I. INTRODUCTION 



One of the main motivations to study high energy 
heavy ion collisions is the creation of a new decon- 
fined phase of strongly interacting matter, the so called 
Quark-Gluon Plasma (QGP) [JI2]. At the Relativistic 
Heavy Ion Collider (RHIC) many experimental observa- 
tions like, e.g., jet quenching and high elliptic flow hint 
to the fact that a strongly coupled QGP (sQGP) might 
have been created 0, i, i, @]. At CERN-SPS energies 
evidence for the creation of a new state of matter has 
been published, e.g., the enhanced K/ir ratio ('horn') and 
the step in the mean transverse mass excitation function 
for pions, kaons and protons [7J. Especially the low en- 
ergy (high /is) program at SPS showed a culmination 
of exciting results. Therefore this energy regime will be 
the subject to further detailed studies at the CERN-SPS, 
BNL-RHIC, JINR-NICA and at the FAIR facility. 

Since the direct detection of free quarks and gluons 
is impossible due to the confining nature of QCD, it is 
important to model the dynamical evolution of heavy 
ion reactions to draw conclusions from the final state 
particle distributions about the interesting early stage 
of the reaction. One approach which aims at the de- 
scription of heavy ion reactions consistently from the 
initial state to the final state is relativistic transport 
theory @, i, EE EH, El El • This microscopic descrip- 
tion has been applied quite successfully to the partonic 



as well as to the hadronic stage of the collision. Un- 
fortunately, most transport approaches are restricted to 
2 — > n scattering processes. Thus, if the particle density 
increases it becomes questionable if a restriction to two- 
particle interaction is still justified. While first attempts 
to include multi-particle interactions have been proposed 
[TTI . E3, EE EE E3|, this extension of transport theory 
is still in its infancy. To explain hadronization and the 
phase transition between the hadronic and the partonic 
phase on a microscopic level is also one of the main open 
issues that still has to be resolved. It is therefore difficult 
to find an appropriate prescription of the phase transition 
in such a microscopic approach. First, however promising 
attempts to solve the microsco pic hadronizati on p roblem 
can be found in the literature [H El EE EE Eli El • 

Hydrodynamics, on the other hand, has been proposed 
many years ago as a tool for the description of the hot 
and dense stage of heavy ion reactions where the mat- 
ter might behave like a locally thermalized ideal fluid 
[IE El, 25]- In this approach it is possible to model phase 
transitions explicitly because one of the major inputs 
to a hydrodynamic calculation is the equation of state 
(EoS). The hydrodynamic description has gained impor- 
tance over the last few years because the high elliptic flow 
values that have been observed at RHIC seem compatible 
with some ideal hydrodynamic predictions [IE EE Ell- 
The initial conditions and freeze-out prescription are the 
boundary conditions for a hydrodynamic calculation and 
therefore a further crucial input. Thus, the hydrody- 
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namic results depend strongly on the initial and final 
state prescription that is applied in the specific calcula- 
tion. 

To get a more consistent picture of the whole dynam- 
ics of heavy ion reactions various so called microscopic 
plus macroscopic (micro+macro) hybrid approaches have 
been launched during the last decade. Most notewor- 
thy in this respect are the pioneering studies related to 
a coupling between UrQMD (Ultra-relativistic Quantum 
Molecular Dynamics) and hydrodynamics (a detailed sys- 
tematic investigation of this coupling procedure can be 
found in the following references [13, H3, [33, US [H, [H, 

Other approaches in the same spirit are, e.g., the NEX- 
SpheRIO approach that uses initial conditions calculated 
in a non-equilibrium model (NEXUS) followed by an 
ideal hydrodynamic evolution [1| [IS S [M S3, El El 
or a hybrid approach by Toneev et al. which uses QGSM 
initial conditions followed by a three-dimensional hy- 
drodynamic evolution [43| . In this way event-by-event 
fluctuations are taken into account and the calculation 
mimics more realistically the experimental case. For the 
freeze-out NEXspheRIO employs a continuous emission 
scenario or a standard Cooper-Frye calculation. Other 
groups, e.g., Teaney et al. Hirano et al. [H, S3], 

Bass/Nonaka [25|, are using smooth Glauber or Color 
Glass Condensate initial conditions followed by a full two- 
or three-dimensional hydrodynamic evolution and calcu- 
late the freeze-out by a subsequent hadronic cascade. 
The separation of chemical and kinetic freeze-out and 
final state interactions like resonance decays and rescat- 
terings arc taken into account. There arc two major con- 
clusions from these previous studies: The treatment of 
the initial state fluctuations and the final decoupling is 
of major importance for a sound interpretation of the 
experimental data. 

Unfortunately, all presently existing micro+macro ap- 
proaches rely on a complete separation of the three main 
ingredients (initial conditions, hydrodynamic evolution, 
transport calculation). Thus, it is impossible to com- 
pare the evolution of the system between hydrodynamics 
and transport simulation directly and from the same ini- 
tial conditions. This may provide essential new insights 
into the role of viscosity and local equilibration. In addi- 
tion, the usual separation of the program code does not 
allow for a dynamical coupling between hydrodynamics 
and transport calculation, which would be desirable to 
consistently solve the freeze-out puzzle [IS S3, S3, H3| • 

To overcome these restrictions, we go forward and 
build a transport approach with an embedded three- 
dimensional ideal relativistic one fluid evolution for the 
hot and dense stage of the reaction. This allows to re- 
duce the parameters for the initial conditions and the 
freeze-out prescription. The aim is to compare calcula- 
tions with different EoS within the same framework. It 
will be possible to extract the effect of changes in the 
EoS, e.g., a phase transition from hadronic matter to the 
QGP, on observables. 



In this paper we describe the specific micro+macro 
hybrid approach that embeds a hydrodynamic phase in 
the UrQMD approach. First we explain the initial con- 
ditions, then introduce the basics of the hydrodynamic 
evolution including the hadron gas EoS and the transport 
calculation and illustrate how the freeze-out is treated. 
In the second part, the sensitivity of the results on the pa- 
rameters are tested and the time evolution of the baryon 
density and the particle numbers are compared. Results 
on particle multiplicities and rapidity as well as trans- 
verse mass spectra are presented in the third part. 

At present we have calculated results imposing a 
hadron gas EoS to provide a baseline calculation to 
disentangle the effects of the different assumptions for 
the underlying dynamics in a transport vs. hydrody- 
namic calculation. The purely hadronic calculations can 
be compared in the broad energy regime from £i a b = 
2 — 16CL4 GeV where a vast amount of experimental data 
from BNL-AGS and CERN-SPS exists and which will 
be explored in more detailed energy scans by the FAIR 
project near GSI and the critRHIC program. Studies 
employing different EoS are delayed to future work to 
concentrate on effects of the underlying dynamics first. 



II. GENERAL ASPECTS 

The modelling of the dynamical evolution of heavy ion 
reactions is essential to gain further insights about the 
properties of the newly produced hot and dense QCD 
matter. Transport theory aims at the description of all 
stages of the collision on the basis of an effective solution 
of the relativistic Boltzmann equation |5l| 

^•S M /i(s v J p")=C f . (1) 

This equation describes the time evolution of the dis- 
tribution functions for particle species i and includes the 
full collision term on the right hand side. The interaction 
with external potentials leads to an additional term on 
the left hand side. The influence of potentials gets small 
at higher energies compared to the energy that is trans- 
ferred by collisions. Therefore, they are dropped in Eqn. 
[Hand are not further discussed here. Usually, the colli- 
sion kernel is truncated on the level of binary collisions 
and 2 —* n processes to keep the calculation numerically 
tractable. This microscopic approach has the advantage 
that it is applicable to non-equilibrium situations and the 
full phase space information is available at all stages of 
the heavy ion reaction. The restriction to binary colli- 
sions assumes large mean free paths of the particles. Be- 
tween interactions the particle trajectories are given by 
straight line trajectories and particles are assumed to be 
in asymptotic states between the collisions (no "memory 
effect"). 

This assumption might not be justified in the hot and 
very dense stage of heavy ion collisions anymore. In this 
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regime the continuum limit in form of relativistic hydro- 
dynamics might fit better to the characteristics of the 
system. The hydrodynamic evolution is governed by the 
energy and momentum conservation laws for given ini- 
tial conditions, i.e. spatial distributions of energy and 
net baryon number densities. The coordinate space is 
divided into small cells in which the distribution func- 
tions correspond to equilibrium distributions (Fermi or 
Bose distribution). In this macroscopic approach the 
propagated quantities are net baryon number and energy 
densities which can be translated into information about 
the temperature and chemical potential via the specific 
equation of state (EoS). Since the evolution is driven by 
pressure gradients and the pressure is determined via the 
EoS, the EoS is the essential ingredient for the hydrody- 
namical evolution. Thus, hydrodynamics is a good tool 
to describe collective behaviour. Ideal hydrodynamics 
applies to systems with small mean free path, otherwise 
viscous effects have to be taken into account [52j . A gen- 
eral advantage of hydrodynamics is the feature to explic- 
itly incorporate phase transitions by changing the EoS. 

However, in the late stage of the heavy ion reaction 
the system gets too dilute to apply ideal fluid dynamics. 
The hadronic rescatterings and decays of resonances have 
to be described, e.g., by using a transport description. 
Overall, there are two crucial points one has to take care 
of when building up a transport+hydrodynamics hybrid 
approach. The first is the initial switch from the mi- 
croscopic to the macroscopic calculation where it has to 
be ensured that the local equilibrium assumption is ful- 
filled. The second one is the so called freeze-out where 
the hydrodynamic fields are mapped to particles that are 
further propagated in a hadronic cascade. The freeze-out 
transition must be placed in a region where both descrip- 
tions are valid at the same time, e.g., the phase transition 
region. In the following the specific implementation de- 
veloped here will be discussed in more detail. 



III. SPECIFIC MICRO+MACRO APPROACH 

A. UrQMD Approach 

For our investigation, the UrQMD model (v2.3) [l©] is 
applied to heavy ion reactions from E} a \, = 2— 160 A GeV. 
This non-equilibrium transport approach constitutes an 
effective solution of the relativistic Boltzmann equation 
(see Eqn. [1]). The underlying degrees of freedom are 
hadrons and strings that are excited in high energetic bi- 
nary collisions. Mean fields can in principle be taken into 
account in this framework, but the model is run in the 
so called cascade mode without inter-particle potentials. 
To omit the potentials is reasonable because the inclu- 
sion of mean fields would not change the results in the 
energy range that we are considering here. Note that this 
is consistent with the calculation of the equation of state 
for the hydrodynamic evolution where no mean field has 
been taken into account as well. 



The projectile and target nuclei are initialised accord- 
ing to a Woods-Saxon profile in coordinate space and 
Fermi momenta are assigned randomly for each nucleon 
in the rest frame of the corresponding nucleus. The 
hadrons are propagated on straight lines until the col- 
lision criterium is fulfilled. If the covariant relative dis- 
tance dtrans between two particles gets smaller than a 
critical distance that is given by the corresponding total 
cross section a collision takes place, 



dtrans < ^0 = J , Otot = ^(v^, type) . (2) 

Each collision process is calculated in the rest frame of 
the binary collision. The reference frame that is used for 
the time ordering of the collisions and later on also for the 
switchings to and from the hydrodynamic phase is the 
equal speed-system of the nucleus-nucleus collision (for 
symmetric systems the equal speed system is identical to 
the center of mass system). 

In UrQMD 55 baryon and 32 meson species, ground 
state particles and all resonances with masses up to 
2.25 GeV, are implemented with their specific properties 
and interaction cross sections. In addition, full particle- 
antiparticle symmetry is applied. Isospin symmetry is 
assumed and only flavour-SU(3) states are taken into ac- 
count. The elementary cross sections are calculated by 
detailed balance or the additive quark model or are fitted 
and parametrized according to the available experimen- 
tal data. For resonance excitations and decays the Breit- 
Wigner formalism, utilizing their vacuum properties is 
employed. 

Towards higher energies, the treatment of sub- hadronic 
degrees of freedom is of major importance. In the present 
model, these degrees of freedom enter via the introduc- 
tion of a formation time for hadrons produced in the frag- 
mentation of strings [H, HH, [55| . String excitation and 
fragmentation is treated according to the Lund model. 
For hard collisions with large momentum transfer (Q > 
1.5 GeV) PYTHIA is used for the calculation. A phase 
transition to a quark-gluon state is not incorporated ex- 
plicitly into the model dynamics. However, a detailed 
analysis of the model in equilibrium yields an effective 
equation of state of Hagedorn type [Eg, H3, [H, [EH, [6(| • 
The UrQMD transport model is successful in describ- 
ing the yields and the pt spectra of various particles in 
proton-proton, proton-nucleus and nucleus-nucleus colli- 
sions [61j . A compilation of results of the actual version 
UrQMD-2.3 compared to experimental data can be found 
in [63- 

Apart from the success of transport simulations to de- 
scribe spectra and yields certain problems remain: 

• Elliptic flow values above SPS energies are too 
small jHHi, 

• HBT radii hint to a very small R /R s ratio [65l.l66l]. 

• Strangeness, especially multi-strange baryons are 
not produced in sufficient amounts [67j |. 
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These observables that are sensitive to the early stage 
of the collision (pressure) or to the approach of ther- 
mal and chemical equilibrium during the collision history 
hint to the fact that a purely hadronic transport model 
may not be sufficient to describe the dynamics of the 
hot and dense stage of heavy ion reactions at higher en- 
ergies [H, [64], [H, lH, [6i| . Therefore, these observations 
exemplify the need to embed a full three-dimensional rel- 
ativistic fluid dynamics description for these stages of the 
reaction. 

For the results that are shown in this paper, the ref- 
erence calculations are always performed employing the 
state-of-the-art UrQMD-2.3 model. 



B. Initial Conditions 

The Ultra-relativistic Quantum Molecular Dynamics 
Model is used to calculate the initial state of a heavy ion 
collision for the hydrodynamical evolution [35| . This is 
necessary to account for the non-equilibrium nature of 
the very early stage of the collision. Event-by-event fluc- 
tuations of the initial state are naturally included by this 
set-up. The coupling between the UrQMD initial state 
and the hydrodynamical evolution takes place when the 
two Lorentz-contracted nuclei have passed through each 
other. The initial time to begin with the hydrodynamical 
evolution is calculated via Eqn. [3] (and is assumed to be 
at least 1 fm/c): 



2R 



2R 



= 2R 



lab 



(3) 



where R is the radius of the nucleus, mpj is the nu- 
cleon mass and -E\ab is the kinetic beam energy. This 
assures that (essentially) all initial baryon-baryon scat- 
terings have proceeded and that the energy deposition 
has taken place. This is the earliest possible transition 
time where thermalization might be achieved 60]. It is 
also convenient from the hydrodynamical point of view 
since at that time the two baryon currents that fly into 
oppposite directions have separated again. 

In general, it is not well-established how and when 
chemical and kinetic equilibrium might have been 
reached in the early stage of the collision. One of the 
problems is, e.g., that the local equilibrium assumption 
might not apply equally well to all parts of the system at 
the same time in the computational frame which corre- 
sponds to the center of mass system of the two colliding 
nuclei. As a consequence, the faster particles have had 
less time in their local rest frame to equilibrate. For 
the bulk part and the high density region at midrapidity 
the difference between the two frames is small. These 
problems are present in all hydrodynamic/macroscopic 
approaches that rely on an equilibrium assumption and 
it is not our attempt to resolve these difficulties in this 
paper. One perspective might be the dynamical coupling 



between the initial transport calculation and the hydro- 
dynamic evolution including source terms on both sides 
of the transition surface. 

To allow for a consistent and numerically stable map- 
ping of the 'point like' particles from UrQMD to the 3- 
dimensional spatial-grid with a cell size of (0.2 fm) 3 , each 
hadron is represented by a Gaussian with a finite width. 
"Pre-formed" hadrons in the process of string fragmen- 
tation are also included in the transformation to the hy- 
drodynamic quantities. I.e. each particle is described by 
a three-dimensional Gaussian distribution of its total en- 
ergy, momentum (in x-, y-, and z-direction) and baryon 
number density. The width of these Gaussians is cho- 
sen to be <j = l fm. A smaller Gaussian widths leads 
to numerical instabilities (e.g., entropy production) in 
the further hydrodynamical evolution, while a broader 
width would smear out the initial fluctuations to a large 
extent. To account for the Lorentz-contraction of the 
nuclei in the longitudinal direction, a Lorentz-gamma- 
factor is included. The resulting distribution function in 
the computational frame (cf), e.g., for the energy density, 
reads: 



e c f(x,y,z) = Ne~ 



i P ) +(y-y P ) +(72(2 
^ 



(4) 



where N = (l/2n) 3 / 2 "( z /(j 3 E c f provides the proper nor- 
malisation, e c f and E c f are the energy density and total 
energy of the particle in the computational frame, while 
(xp,y p , z p ) is the position vector of the particle. Sum- 
ming over all single particle distribution functions leads 
to distributions of energy, momentum and baryon num- 
ber densities in each cell. 

To allow for calculations at finite impact parameter the 
spectators - nucleons that have not interacted at all be- 
fore the start time of the hydrodynamic evolution t star t 
- are propagated separately from the hydrodynamic evo- 
lution. The spectators are propagated on straight line 
trajectories in the usual cascade mode until the end of 
the hydrodynamic phase has been reached. 

Instead of smearing out the initial distributions by de- 
scribing the point like hadrons as Gaussian distributions, 
one could also obtain a smooth distribution by averaging 
over a large sample of UrQMD events. Our procedure of 
creating a new initial state for each event is motivated 
by the fact, that the experimental results all relate to 
observed (averaged) final, and not initial, states. Thus, 
event-by-event fluctuations of the initial state can be ob- 
served, (e.g., in V2 fluctuations) and have therefore been 
taken into account properly (for discussion of the impor- 
tance of these fluctuations see, e.g., [36l. l70j). 

As an example, Figs. [T] and [2] show the energy and 
baryon number densities obtained in one single central 
(6 = fm) Pb+Pb collision at £ lab = 40,4 GcV after 
the initialisation of the hydrodynamic fields. The start- 
ing time is in this case t s t a rt — 2.83 fm and the den- 
sities in the figures correspond to the same time. All 
quantities are given in the local rest frame. The max- 
imum values reach 6 GeV/fm 3 for the energy density 
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2.83 fm Energy density distribution 




FIG. 1: (Color online) Initial energy density distribution in 
the reaction plane (x — z plane) of one central (6 = fm) 
Pb+Pb collision at E\ab = 40^4 GeV. z corresponds to the 
beam direction and x to the in-plane axis (direction of the 
impact parameter) of the collision. 

2.83 fm Net baryon number density distribution 



ergy momentum tensor which is also known as the Lan- 
dau frame. This frame coincides in ideal hydrodynamics 
with the Eckart frame which is defined as the local rest 
frame of the baryon number current. The iterative cal- 
culation of the cell velocity succeeds if those two frames 
are close enough to each other. By this transformation 
the system is forced to local equilibrium. 

C. Hydrodynamic Evolution 

Ideal relativistic one fluid dynamics is based on the 
conservation of energy, momentum and the net baryon 
number current. For the hydrodynamical evolution local 
equilibrium is assumed and zero viscosity which corre- 
sponds to zero mean free path. The two conservation 
equations that govern the evolution are [12, [Zll 

9 M T^ = and d^N* = 0, (5) 

where is the energy-momentum tensor and is the 
baryon current. For an ideal fluid the energy- momentum 
tensor and the net baryon number current take the simple 
form 




FIG. 2: (Color online) Initial net baryon number density dis- 
tribution in the reaction plane (x — z plane) of one central 
(6 = fm) Pb+Pb collision at Ei^ = 40^4 GeV. z corresponds 
to the beam direction and x to the in-plane axis (direction of 
the impact parameter) of the collision. 

and around 8 times the nuclear ground state density 
for the baryon number density. The distributions are 
quite smooth which is necessary to provide possible ini- 
tial conditions for the hydrodynamic evolution. One can 
see some peaks that correspond to local maxima of the 
distributions ("hot spots"). It is further clear that the 
single event distributions are not symmetric, neither in 
the transverse nor the longitudinal direction. 

The remaining question is if the system is thermalizcd 
enough to assure that the local equilibrium assumption 
of ideal hydrodynamics is fulfilled. In our case, the hy- 
drodynamic code transforms all the given quantities from 
the computational frame to the local rest frame of the en- 



T" v = (ejrf + P) u v - P g^ and = p M w M (6) 

where ei r f,-P and p\ T { are the local rest frame energy 
density, pressure and net baryon density, respectively, 
u'' = j(l,v) is the four velocity of the cell and g^ v — 
diag(+, — , — , — ) is the metric tensor. The local rest 
frame is defined as the frame where T^ v has diagonal 
form, (i.e. all off-diagonal elements vanish). The four- 
velocity of the cells is calculated via the transformation 
into the local rest frame. 

The equations of motion are solved in the following 
form by employing computational frame quantities e c f ,p l 
and p c f for the energy, momentum and net baryon num- 
ber densities. 

d t e c{ + V -{ecfv) = -V ■ (Pv) (7) 
d t p + V ■ {pv) = -VP (8) 
dtPot +V- (p ci v) = (9) 

In our case, the full (3+1) dimensional hydrodynamic 
evolution is performed using the SHASTA algorithm 
[72I [lH . The partial differential equations are solved on 
a three-dimensional spatial Eulerian grid with fixed posi- 
tion and size in the computational frame. The standard 
size of the grid is 200 cells in each direction while the cell 
size has been chosen to be (da:) 3 = (0.2 fm) 3 which leads 
to timesteps of dt = 0.08 fm. Depending on the beam 
energy, the cell sizes may require adjustment to assure a 
stable solution of the differential equation. 

The equation of state is needed as an additional input 
to calculate the pressure, temperature and chemical po- 
tential corresponding to the energy and the baryon num- 
ber densities. Since the evolution of the system is driven 
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by pressure gradients the EoS has the most important 
influence on the evolution. 



D. Equation of State 

To solve the hydrodynamical equations, the EoS, the 
pressure as a function of energy and net-baryon number 
density, is needed as an input. Since the actual EoS of hot 
and dense QCD matter is still not precisely known, it may 
seem disadvantageous to have this additional uncertainty 
in the model. On the contrary it may prove to be an im- 
portant trait of the model to be able to study changes 
on the dynamics of the bulk matter when changing the 
EoS thus finding observables for a phase transition in hot 
QCD matter. For recent discussions of different EoS and 
how to obtain EoS from lattice calculations, the reader 
is referred to 0, M, ~Zfj. 

The EoS used in the present calculations is a grand 
canonical description of a free, non interacting gas of 
hadrons. We will refer to it as the Hadron Gas (HG). 
It follows from the hadronic chiral model presented in 
[77l . l78j . The chiral hadronic SU (3) Lagrangian incorpo- 
rates the complete set of baryons from the lowest flavour- 
SU(3) octet, as well as the entire multiplets of scalar, 
pseudo-scalar, vector and axial- vector mesons [7""] . In 
mean-field approximation, the expectation values of the 
scalar fields relevant for symmetric nuclear matter corre- 
spond to the non-strange and strange chiral quark con- 
densates, namely the a and its ss counterpart £, respec- 
tively, and further the u> and 4> vector meson fields. An- 
other scalar iso-scalar field, the dilaton x, is introduced 
to model the QCD scale anomaly. However, if \ does 
not couple strongly to baryonic degrees of freedom then 
it remains essentially "frozen" below the chiral transition 
[7""] . Consequently, we focus here on the role of the quark 
condensates. 

Interactions between baryons and scalar (BM) or vec- 
tor (BV) mesons, respectively, are introduced as 

£bm = - A {giao + 9tcO i>< , (10) 

i 

£bv = - ^2 tpi (g^ou + gi<t>io4>°) ipi , (11) 

i 

Here, i sums over the baryon octet (N, A, S, S). A term 
£ vcc with mass terms and quartic self-interaction of the 
vector mesons is also added: 



£vec = ^a^x 2 ^ 2 + ^a<f>X 2 2 + 34 (^ 4 + 2( ? !)4 ) • 



.2. .2 1 



The scalar self-interactions are 

Co = -^ 0X 2 (a 2 +C 2 ) + fcl^ 2 + C 2 ) 2 + fc2(y+C 4 ) 

+ k 3X a\ ~ fc 4X 4 - \ X A In 4 + U 4 ln 4t- • (I 2 ) 
4 Xo 3 °oCo 



Interactions between the scalar mesons induce the spon- 
taneous breaking of chiral symmetry (first line) and the 
scale breaking via the dilaton field x (last two terms). 

Non-zero current quark masses break chiral symmetry 
explicitly in QCD. In the effective Lagrangian this corre- 
sponds to terms such as 



r x 

i-SB — 2 

Xo 



mlfn<7 + (V2m 2 K f K - -j=™l M( 



(13) 



According to £bm f"U|) . the effective masses of the 
baryons, m*(a, C) = gi a & + g^ C , are generated through 
their coupling to the chiral condensates, which attain 
non-zero vacuum expectation values due to their self- 
interactions ("""J in Co (|12|) . The effective masses of the 
mesons are obtained as the second derivatives of the 
mesonic potential Vmosoii = — Co — C vcc — Csb about its 
minimum. 

All parameters of the chiral model discussed so far 
are fixed by either symmetry relations, hadronic vacuum 
observables or nuclear matter saturation properties (for 
details see |77[). In addition, the model also provides a 
satisfactory description of realistic (finite-size and isospin 
asymmetric) nuclei and of neutron stars [zU H3] • 

If the baryonic degrees of freedom are restricted to the 
members of the lowest lying octet, the model exhibits 
a smooth decrease of the chiral condensates (crossover) 
for both high T and high /j,b [ZZl [Sll- However, addi- 
tional baryonic degrees of freedom changes this into a 
first-order phase transition in certain reg imes of the T- 
Hb plane, depending on the couplings 35l[78l. l8ll. [82* 



In what follows, the meson fields are replaced by 
their (classical) expectation values, which corresponds to 
neglecting quantum and thermal fluctuations. Fermions 
have to be integrated out to one-loop. The grand 
canonical potential can then be written as 



Q./V — —C vcc — Co — Csb — V v . 



(14) 



(2tt) 3 

li 
(2tt) 3 



where 7b, 7a/ denote 
spin-isospin degeneracy 



d A k 



ln(l + e-TK«-ft*: 



ml 



the baryonic 
factors and 



and 
E 



B.M 



mesonic 
(fc) = 



k 2 



l B,M 



are the corresponding single particle en- 



ergies. The effective baryon-chemical potentials are /i* = 
fJ*i-giuU-9i<j><f>, with fii = (7i % q -n\)ii q + (n l s ~7i\)ii s . Here 
Mg = Mb/3 is the quark, and /i s the strange quark, chem- 
ical potential. The potentials of the mesons are given by 
the sum of the corresponding quark and anti-quark chem- 
ical potentials. The vacuum energy V vac (the potential 
at pB = T = 0) has been subtracted. 

By extremizing tt/V one obtains self-consistent gap 
equations for the meson fields. Here, globally non- 
strange matter is considered and /is for any given T 
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and fiB is adjusted to obtain a vanishing net strangeness. 

Once the grand canonic potential is known as a 
function of T and /is, all other thermodynamic quanti- 
ties are derived straightforward. In its minima the grand 
canonic potential corresponds to —p, the pressure. The 
entropy density s, number density n and energy density 
e then follow from the Euler relation, 



sT 



turn 



(15) 



where the sum runs over all included hadron species. 

Setting all hadron masses and chemical potentials 
to their vacuum values, and adding all reliably known 
heavy resonance states - with masses up to 2 GeV [H| 
- as free particles into (fF4"j) . yields the above mentioned 
Hadron Gas EoS [85j]. Hence, the hadronic degrees of 
freedom included in this EoS are consistent with the 
active degrees of freedom in the UrQMD model. This 
enables us to directly compare the dynamics of the 
hydrodynamic model with the transport simulation. 

Using the chiral model and adding additional baryonic 
degrees of freedom as well as adjusting their scalar and 
vector coupling, an EoS with a phase structure including 
a first order phase transition and even a critical endpoint 
at finite [Ib can be obtained [78j]. This chiral EoS has 
already been successfully applied to a hydrodynamic 
calculation (35[. Here the essentially different equation 
of state leads to distinguishable different results on the 
properties of bulk matter. 

To emphasize these differences even more, a MIT-Bag 
model EoS matched to an interacting hadron gas [7^ |. 
generating a phase structure with a broad first order 
phase transition at all /is, can also be applied in our 
model. 

Consequently these results will be compared to our 
present calculations, constituting observables for a phase 
transition in hot and dense QCD matter, even suggesting 
the order of this phase transition. Comparisons between 
the different EoS (i.e. free Hadron Gas, chiral EoS, and 
Bag model EoS) will be presented in a follow-up paper. 



(16) 



where f(x,p) are the boosted Fermi or Bose distri- 
butions corresponding to the respective particle species. 
Since we are dealing with an isochronous freeze-out, the 
normal vector on the hypersurface is da^ — (d 3 x,0). 

Let us note that it is of utmost importance to con- 
sider the same degrees of freedom on both sides of the 
hypersurface because otherwise energy and momentum 
conservation is violated. In our case, this is assured by 
the inclusion of the same particle species in the equa- 
tion of state for the hydrodynamic calculation as in the 
transport calculation. In principle, it might also happen 
that particles are moving back into the hydrodynamic 
phase, however, the explosive character of heavy ion re- 
actions, i.e. the rapid expansion flow suppresses the back- 
streaming effect. Therefore, this effect is negligible in our 
situation [871 ] . 

The assumption of an isochronous freeze-out leads to 
fluctuations of the temperature and baryo-chemical po- 
tential distributions and not to single values for some 
of the thermodynamic quantities on the hypersurface. 
The quality of this approach and that the two (hydro- 
dynamics and transport) prescriptions are valid through 
the applied range of switching temperatures is shown in 
the next section, where parameter tests and the time evo- 
lutions of the particle yields are explored in detail. 
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E. Freeze-out 



FIG. 3: (Color online) Distribution of the energy in the cells 
at freeze-out at _Ei a b = 40^4 GeV. 



Presently, the hydrodynamic evolution is stopped, if 
the energy density drops below five times the nuclear 
ground state energy density (i.e. ~ 730 MeV/fm 3 ) in all 
cells. This criterium corresponds to a T-^s-configuration 
where the phase transition is expected (see dotted line in 
Fig. Uj), i.e. a region where the hydrodynamic and the 
transport description are valid at the same time. The 
hydrodynamic fields are mapped to hadrons according 
to the Cooper-Frye equation [86| 



Fig. [3] shows the distribution of energy in the cells 
at freeze-out with respect to temperature and baryo- 
chemical potential at = 40A GeV. We present here 
the energy distribution and not the number of cells be- 
cause it is more interesting to know where the energy of 
the system sits than considering the many almost empty 
cells that do essentially not contribute to particle pro- 
duction. From Fig. [3] one obtains the mean values of the 
distributions that are in line with results from statistical 
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model fits. The mean values of, e.g., the temperature can 
be calculated as 



(T) 



E,,j,fc T(i,j,k)p B (i,j,k) 



(17) 



where i,j, k are the cell indices and the sum runs over 
all cells of one event. The net baryon number density ps 
has been used as a weighting factor. 
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FIG. 4: (Color online) Mean values of temperatures and 
baryo-chemical potentials at freeze-out for different beam en- 
ergies are depicted as red circles (starting in the lower right 
corner at Ei ah = 2A GeV and going through 4,6,8,11,20,30,40 
and 80 to _Ei ab = 160 A GeV in the upper left). The error 
bars indicate the width of the distribution. The dotted line 
depicts the line of constant energy density (e = 5 ■ eo) that 
corresponds to our freeze-out criterium. For comparison the 
freeze-out line calculated by Cleymans et al. [88J,|89| (full line) 
and results from Dumitru et al. 01 (green open squares) are 
shown. 

To illuminate this finding in more detail, Fig. 2] shows 
the mean values of the temperature and baryo-chemical 
potential distributions at different energies in the T—ps- 
plane for central (6 = fm) Au+Au/Pb+Pb collisions. 
Also, the widths of the distributions are depicted as "er- 
ror" bars. Fig. |4]shows that the present freeze-out distri- 
butions are similar to the parametrized curve for chemical 
freeze-out as calculated by Cleymans et al. from statis- 
tical model fits to final particle multiplicities. The cal- 
culation by Dumitru et al. shows mean values as well as 
widths of temperature and baryo-chemical potential dis- 
tributions that have been obtained by statistical model 
fits to final particle yields employing the assumptions of 



an inhomogeneous freeze-out hypersurface. This calcula- 
tion also leads to similar values as our calculation. 

The effect that the mean temperature at the transi- 
tion to the transport prescription saturates or even drops 
down a little at higher beam energies is related to the 
rapidity distribution of the temperature in the hydrody- 
namic cells at freeze-out which is shown in Fig. Ofor three 
different beam energies. At low beam energies the midra- 
pidity region coincides with the hottest region at freeze- 
out. At higher SPS energies the situation changes. The 
hottest cells are at high rapidities while the midrapidity 
region has already cooled down well below the tempera- 
ture of 170 MeV. This problem might be resolved by a 
different freeze-out prescription on another hypersurface 
(e.g., isotherm, iso-e) and is subject to future investiga- 
tions. The best solution will be the dynamical coupling 
between hydrodynamics and transport which allows also 
for back-streaming contributions. 
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FIG. 5: (Color online) Rapidity profile of the freeze-out tem- 
peratures in the spatial plane with x — y = fm for cen- 
tral Au+Au/Pb+Pb collisions at three different beam ener- 
gies (Slab = 6, 40 and 160A GeV). 

In the following the practical implementation will be 
explained in more detail. The implementation is based 
on a Monte Carlo sampling of Eqn. [TJ] and follows the 
general steps: 

1. The particle numbers Ni are calculated according 
to the following formula, 



Ni 



7 ' Vceii = / d 3 pfi(x,p) ■ 7 • V C( 



(18) 



where the index i runs over the different particle 
species like, e.g., tt, p, p or A. 7 is the boost fac- 
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tor between the computational frame and the cell. 
V ce u is the volume of the cell in the computational 
frame and n is the particle number density. All cells 
with temperatures that are lower than 3 MeV are 
discarded from the following procedure because of 
numerical reasons. The local rest frame equilibrium 
distribution function is denoted by fi(x,p). To sim- 
plify the calculation, a relativistic Boltzmann dis- 
tribution is used for all particles, except pions. It 
has been checked, that the Boltzmann approxima- 
tion is sufficient to describe all particle species it 
is applied to. For the Boltzmann distribution the 
momentum integration leads to the following result 
for the particle number density 

*=-|^«p(£Mt) (19) 

where g is the degeneracy factor for the respective 
particle species, m is the mass of the particle to 
be produced, T the temperature of the cell and K2 
is the modified Bcsscl function. The chemical po- 
tential [i includes the baryo-chemical potential and 
the strangeness chemical potential in the following 
way 

/1 = B ■ [i B + S ■ fi s (20) 

where S is the quantum number for strangeness and 
B is the baryon number. 

For pions the Bose distribution has to be taken into 
account because the pion mass is on the order of the 
temperature of the system. In this case, the mo- 
mentum integration involves an infinite sum over 
modified Bessel functions 

To calculate the number of particles in the compu- 
tational frame the particle number density has to 
be multiplied with the Lorentz-stretched volume of 
the cell (Vceu = (0.2) 3 fm 3 ). 

2. The average total number of particles in the cell, 
(A r ), is the sum over all particle numbers A,; = 

<A)=]Ta,. (22) 

i 

3. The total number of particles emitted from a cell, 
Ni, is obtained from a Poisson distribution accord- 

ing to P(N) = ^e-W. 



In the limit of small mean values, the Poisson dis- 
tribution becomes P(l) ~ (A). Thus it can be 
decided by one random number between and 1 
if a particle is produced in the respective cell. If 
the random number is smaller than (N) one parti- 
cle is produced and there is no particle production 
otherwise. The full Poisson distribution is used, if 
the particle number (A) is larger than (0.01). This 
assures an accuracy better than 1 %. 

4. The particle type is chosen according to the prob- 
abilities Ni/{N}. 

5. The 73 component of the isospin is distributed ran- 
domly because UrQMD assumes full isospin sym- 
metry. To conserve the overall charge of the system 
and the initial isospin-asymmetry the probability to 
generate the isospin component that leads to the 
right value of the charge that should be obtained 
in the end is favoured. The other isospin compo- 
nents are exponentially suppressed. The power of 
the exponential is proportional to the difference of 
the total charge generated by this produced particle 
and the required value. 

6. The 4-momenta of the particles are generated ac- 
cording to the Cooper- Frye equation (see Eqn. 1 16[) . 
For baryons and strange mesons the chemical po- 
tentials for baryon number and strangeness are 
taken into account. 

7. The particle vector information is transferred back 
into the UrQMD model. The subsequent hadronic 
cascade calculation incorporates important final 
state effects as, e.g., rescatterings of the particles 
and resonance decays. 

The above mentioned steps are pursued on random 
cells until the initial net baryon number is reached. 
Strangeness and charge are also conserved in each event 
separately, energy conservation is fulfilled for the mean 
values averaged over several events. Aiming at a real- 
istic description of heavy ion reactions we perform the 
freeze-out for each event separately and do not average 
over many freeze-outs for one hydrodynamical evolution. 

IV. PARAMETER TESTS 

A. Start-time and Freeze-Out Criterium 

In this section we investigate the dependences of ob- 
servables on parameters of the implementation. Two im- 
portant parameters have to be determined. The first 
one is the starting time t s t a rt which defines the initial 
switch from UrQMD to the hydrodynamic evolution. 
The second parameter is the freeze-out criterium which 
is parametrized as an energy density criterium. While 
varying one parameter we have fixed the other one to the 
default value (1 • t sta rt or 5 ■ eq). 
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FIG. 6: (Color online) Pion and kaon multiplicities (upper 
panel) and mean transverse mass of pions and kaons at midra- 
pidity (\y\ < 0.5) (lower panel) for four different starting times 
for central (6 < 3.4 fm) Au+Au/Pb+Pb collisions. The open 
symbols depict the result at _Ei a b = 11 A GeV while the filled 
symbols are the results at Z?i a b = 40A GeV. 



FIG. 7: (Color online) Pion and kaon multiplicities (upper 
panel) and mean transverse mass of pions and kaons at midra- 
pidity Qy\ < 0.5) (lower panel) for four different freeze-out 
criteria for central (b < 3.4 fm) Au+Au/Pb+Pb collisions. 
The open symbols depict the result at £i a b = 11 A GeV while 
the filled symbols are the results at E\ah = 40^4 GeV. 



Fig. [5] (top) shows calculations of the total, i.e. pion 
and kaon (47r) multiplicities, for four different starting 
times at two beam energies. The open symbols depict 
always the result at the highest AGS energy (£) a b = 
11 A GeV) and the filled symbols are the results at the 
SPS energy (-Ei a b = 40A GeV). The starting time is var- 
ied in factors of the default value that has been calcu- 
lated via Eqn. [3] Displayed are results from halved to 
doubled initial time. One observes a higher pion pro- 
duction for earlier starting times compared to the pion 
production in the standard setup (1 t s tart)- This may 



be explained by the fact that the system is forced more 
strongly to equilibrium and the cascade evolution lasts 
longer. If the hydrodynamic evolution is started at later 
times (1.5 or 2 t start) the resulting pion multiplicities are 
not affected anymore. The kaon yield is essentially not 
sensitive to the switching time. To summarize, varying 
the starting time by a factor of 4 results in a change in the 
pion and kaon production of less than ±10% compared 
to the pion and kaon production in the default configu- 
ration (1 tstart)- In Fig. [S] (bottom) the mean transverse 
mass of pions and kaons at midrapidity is shown. The 
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mean transverse mass values are calculated for the same 
four different starting times at the two exemplary beam 
energies as before. Here the results do also not change 
more than ±15% for a large spectrum of starting times. 
Therefore, our choice of the starting time as the geomet- 
rical criterium when the nuclei have passed through each 
other is sensible and stable. It is the earliest possible 
time where thermalization may have been achieved and 
the baryon currents have disconnected. 

Fig. [7] (top) is the equivalent picture to Fig. [6] (top), 
however displaying the dependence of the total pion and 
kaon multiplicities on the freeze-out criterium. The de- 
fault value for the transition energy density is 5eo while 
we have varied it from (4 — 10) Cq. The higher the freeze- 
out energy density the earlier the hydrodynamic evolu- 
tion is stopped because the cells reach this critical energy 
density value earlier. As a consequence the kaon yields 
raises with an increase of the energy density criterium, 
while the pions remain virtually unchanged for all inves- 
tigated transition criteria. Fig. [7] (bottom) shows the 
results for the pions and kaons mean transverse mass for 
two beam energies and different freeze-out criteria. Again 
one observes only a very weak dependence on the freeze- 
out criterion. Furthermore, the mean transverse mass 
values for the two different meson species are very sim- 
ilar since they acquire the same transverse flow. These 
findings confirm that our choice for the freeze-out cri- 
terium as 5eo is robust. 
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FIG. 8: (Color online) Beam energy dependence of the start- 
ing time tatart (blue full line), the averaged time for the hy- 
drodynamic evolution (thydro) (red dotted line) and the mean 
duration of the hadronic stage (thadronic) (green dashed line) 
of central (6 < 3.4 fm) Au+Au/Pb+Pb collisions. 



B. Timescales 

In this Section the time scales that are important will 
be explored. Let us start with a table that summa- 
rizes the mean durations of the hydrodynamic and the 
hadronic phase of the collision for different starting times 
and freeze-out parameters at = 40A GeV. 



% ' tstart 


x ■ eo 


{thydro} [fm] 


(thadronic) [fm] 


l 


4 


7.68 


15.63 


l 


5 


7.72 


16.07 


l 


6 


6.84 


16.49 


l 


10 


4.60 


17.29 


0.5 


5 


7.03 


14.59 


1.5 


5 


6.22 


17.50 


2 


5 


4.91 


17.61 



TABLE I: This table contains the mean durations of the hy- 
drodynamic evolution and the hadronic calculation afterwards 
for different starting times and freeze-out criteria for central 
(b < 3.4 fm) Pb+Pb collisions at E^ b = 40A GeV. 



The duration of the hydrodynamic evolution is a well- 
defined period for each event because of the isochronous 
freeze-out. The average is therefore an average over 100 



events. The hadronic stage starts when the hydrody- 
namic evolution is stopped and it ends when the particles 
have undergone their last interaction. An interaction can 
be an inelastic or an elastic collision or a decay. 

The averaged time duration of the stages of the reac- 
tion (given in Table[T| reflect the expectations. The lower 
the freeze-out energy density the later the hydrodynamic 
freeze-out proceeds and therefore the hydrodynamic evo- 
lution lasts longer while the hadronic stage is shortened. 
The later the hydrodynamic evolution starts, the bigger 
tstart is, the shorter the hydrodynamic evolution lasts. 
The hadronic phase does not show a clear trend. To first 
approximation the final UrQMD stage lasts for 16.5 ± 2 
fm independent of the parameters. 

Fig. [8] shows the beam energy dependences of the 
timescales for the chosen values of t s t ar t and the freeze- 
out energy density, from E\ a \> = 2 — 160^4 GeV. The 
starting time decreases as a function of beam energy from 
more than 10 fm at low energies to less than 2 fm at 
higher SPS energies. The mean duration of the hydrody- 
namic as well as the hadronic phase of the reaction grow 
with raising energy. 



C. Time Evolutions 

In the following Section we investigate the time evo- 
lution of different quantities and compare the results of 



12 



§5 

CD 

^4 



"1 1 1 1 1 r 

UrQMD 

UrQMD+Hydro(HG) 

center of the grid 

x=y=z=0 f m 




t(fm) 



FIG. 9: (Color online) Time evolution of the net baryon num- 
ber density in the local rest frame for central (6 = fm) 
Pb+Pb collisions at £a ab = 40 A GeV. 



the hybrid model calculation to the UrQMD calculation 
without an hydrodynamic stage. Since the net baryon 
density is directly propagated on the hydrodynamic grid, 
it serves as a good example to compare to the default 
UrQMD calculation. In the microscopic approach the lo- 
cal rest frame density is calculated as the zero component 
of the net baryon number current in the frame where the 
corresponding local velocity vanishes [9lj . 

Fig. [9]shows the time evolution of the net baryon num- 
ber density at the center of a central Pb+Pb collision 
at Si a b = 40,4 GeV. The blue full line depicts the de- 
fault UrQMD calculation while the red full line depicts 
the result of the hybrid model calculation. The result 
has some spikes because here we compare single events. 
There are two important observations: The absolute val- 
ues of the net baryon number densities are very similar in 
both cases and there are no obvious discontinuities at the 
switching points to and from the hydrodynamic model 
calculation. The smoothness of the curve confirms our 
choice of parameters. 

In Fig. [TO] (top and middle) the time evolution of the 
particle yields in the two different models for the most 
central Pb+Pb collisions at E\ a \> = AOA GeV is com- 
pared. The multiplicities at different timesteps are ex- 
tracted from the hydrodynamic evolution by converting 
the number densities to particle numbers via the freeze- 
out procedure (see Section IIII E|) . Fig. [10] (top) depicts 
the total particle multiplicity (red circles and full line) 
and the midrapidity multiplicity (blue squares and full 
line). The full lines indicate the default UrQMD calcu- 
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FIG. 10: (Color online) Time evolution of the total particle 
number and the midrapidity (\y\ < 0.5) yield (upper panel), 
of the total number of pions and nucleons (middle panel) and 
of the conserved quantities (lower panel) for central (6 = fm) 
Pb+Pb collisions at £a ab = 40A GeV. Results of the hybrid 
model calculation UrQMD+Hydro (HG) are depicted with 
symbols while UrQMD-2.3 results are represented by lines. 
The total energy of the system (red circles and line) has been 
divided by eight for visibility reasons. The other conserved 
quantum numbers are net baryon number (orange triangles 
and line), the overall charge (blue squares and line) and the 
strangeness (green diamonds and line). The total strangeness 
(black dots and line) is given by the sum of s- and s-quarks. 



lation, while the symbols show the results of the hybrid 
model. The multiplicities increase rapidly in the initial 3 
fm/c and then decrease a little, followed by a slower con- 
stant raise until the final multiplicity is reached. This 
qualitative behaviour is very similar in both approaches. 
The decrease of the multiplicity can be associated with 
the thermalization because absorption and production 
processes are on the same order. Note again that there 
are no discontinuities at the switching times in the hybrid 
calculation. 

Next, we explore the time evolution for two particle 
species in more detail. In Fig [TO] (middle), the pions (red 
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circles and full line) represent the newly produced par- 
ticles while the nucleons (blue squares and full line) are 
already there in the beginning as the constituents of the 
two incoming nuclei. The qualitative behaviour of the 
temporal evolution of the pion yield is similar to that 
discussed above for the total multiplicity. The decrease 
at the starting time of the hydrodynamic evolution t ~ 3 
fm/c is much stronger than in the model without hydro- 
dynamic phase, because of the instant thermalization at 
the transition time. The default UrQMD transport cal- 
culation results in a similar, but much smoother, shape 
of the curve. This similarity hints to the fact that the mi- 
croscopic calculation also equilibrates the hot and dense 
matter to a rather large degree. The number of nucle- 
ons decreases in the beginning due to the production of 
resonances and string excitations. At the thermalization 
the minimum is reached and the number of nucleons in- 
creases slowly until the final value is reached. In this 
case, not only the qualitative behaviour is independent 
of the underlying dynamics but also the absolute values 
are very close to each other. 

To check quantum number and energy conservation 
the temporal evolution of all the important quantities in 
both approaches for central Pb+Pb collisions at E\ a b = 
40 A GeV is depicted in Fig. [TO] (bottom). The de- 
fault UrQMD calculations are indicated by full lines while 
the hybrid model calculations with integrated hydrody- 
namic evolution are depicted by symbols. The net baryon 
number (orange triangles and full line), the charge (blue 
squares and full line) and the net strangeness (green di- 
amonds and full line) are exactly conserved in both ap- 
proaches. The total energy (red circles and full line) is 
only on average over several events conserved in the hy- 
brid model calculation due to the freeze-out prescription. 
But the fluctuations are on a 5% level. Note however, 
that the total strangeness in the system (s + s— quarks) 
is very different in both approaches. In the default trans- 
port calculation (black line) the total strangeness in- 
creases in the early stage of the collision and remains con- 
stant. This is contrasted by the hybrid calculation (dots). 
Due to the local thermal equilibration and the thermal 
production of strange particles in the hybrid calculation 
the yield of strange quarks jumps to a higher value at 
the switching time (t start)- The total strangeness then 
decreases as the system cools down, but the final value 
remains 50% higher than in the default transport calcu- 
lation. 



D. Final State Interactions 

The last step is the analysis of the freeze-out process, 
i.e. how much hadronic interactions happen after the 
hydrodynamic evolution. For this purpose, we have cal- 
culated the number of collisions during the hadronic cas- 
cade calculation in dependence of time and y/s of the 
elementary collisions. The corresponding distributions 
are shown in Figs. [TTIandfT^l respectively, for three dif- 
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FIG. 11: (Color online) Temporal distribution of binary col- 
lisions in the hadronic cascade calculation after the hydro- 
dynamic evolution. The upper plot depicts the result at 
Biab = UA GeV, the middle plot at Ei ab = 40,4 GeV and the 
lowest plot at the highest SPS energy (Eub = 160A GeV) for 
central (6 < 3.4 fm) Au+Au/Pb+Pb collisions. The full line 
refers to baryon-baryon collisions (B+B), the dotted line to 
baryon-meson collisions (B+M) and the dashed line to meson- 
meson collisions (M+M). The grey shaded area depicts the 
averaged duration of the hydrodynamic evolution. 



ferent beam energies (11, 40 and 160A GeV, from top to 
bottom ). 

Fig. [Tl] shows the collision rates for meson-meson, 
meson-baryon and baryon-baryon interactions. The grey 
area indicates the average time span of the hydrodynamic 
phase. One observes that substantial collision rates are 
present directly after the transition to UrQMD. The col- 
lision rates stay high for 5 fm/c, and after 30-40 fm/c the 
system is completely frozen out. Only some resonance de- 
cays proceed for longer. According to the composition of 
the system the baryon-baryon and baryon-meson interac- 
tions dominate the lower beam energy result, while at the 
highest SPS energy the meson-meson together with the 
meson-baryon interactions are the most abundant type 
of collision indicating the transition from baryon domi- 
nated to meson dominated systems. Note that the over- 
lap of the hadronic interaction phase with the hydrody- 
namic evolution results from the fact that the duration 
of the different stages fluctuates in the present approach. 
Shown here is the average duration of the hydrodynamic 
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FIG. 12: (Color online) Distribution of the y/s values for the 
binary collisions in the hadronic cascade calculation after the 
hydrodynamic evolution. The upper plot depicts the result at 
Slab = 11 A GeV, the middle plot at £a ab = 40 A GeV and the 
lowest plot at the highest SPS energy (2?i a b = 160A GeV) for 
central (b < 3.4 fm) Au+Au/Pb+Pb collisions. The full line 
refers to baryon-baryon collisions (B+B), the dotted line to 
baryon-meson collisions (B+M) and the dashed line to meson- 
meson collisions (M+M). 



evolution. 

Fig. [T2] shows the ^fs distribution of the elementary 
collision in the freeze-out process. One nicely observes all 
the resonance peaks in the corresponding channels. This 
figure suggests that the most abundant meson-baryon 
collisions are excitations of the A resonance (i.e. ttN 
interactions) since there is a sharp peak at the A mass 
(to a = 1-232 GeV). For meson- meson reactions the p and 
the u> peaks are clearly present. This result indicates that 
there is still resonance regeneration even at this late stage 
of the system's evolution. 



V. RESULTS 

A. Multiplicities and Particle Spectra 

We start with a comparison of the multiplicities and 
particle spectra in the two frameworks. Calculations 
with the embedded hydrodynamic evolution employing a 



hadron gas equation of state for the high density stage of 
the collisions are compared to the reference results of the 
default transport calculation (UrQMD-2.3). Since both 
calculations use the same initial conditions and freeze- 
out prescription it allows to extract, which observables 
are sensitive to the change in the underlying dynamics, 
thus allowing to explore the effect of local equilibration, 
of viscosities and heat conductivity. The hybrid model 
calculation is in the following always depicted as full lines 
while the default UrQMD-2.3 calculations are depicted as 
dotted lines. Note that we do not tune any parameters 
for different energies or centralities for the hybrid model 
calculation. The starting time for the hydrodynamic ex- 
pansion is always calculated using Eqn. [3] and the fixed 
energy density criterium (5eo) for the freeze-out (as ex- 
plained in Section [ill E[) is always employed. 

In Fig. Q2] the excitation functions of the total multi- 
plicities are shown for central Au+Au/Pb+Pb collisions 
for i?i a b = 2 — 160A GeV. The present hybrid approach 
simulations have been restricted to this energy range be- 
cause for calculations at higher energies some numerical 
subtleties have to be resolved, e.g. a dynamical grid size 
for the hydrodynamical evolution. Compared to the de- 
fault simulation, the pion and proton multiplicities are 
decreased over the whole energy range in the hybrid 
model calculation due to the conservation of entropy in 
the ideal hydrodynamic evolution. The non-equilibrium 
transport calculation produces entropy and therefore the 
yields of nonstrange particles are higher. The produc- 
tion of strange particles however, is enhanced due to the 
establishment of full local thermal equilibrium in the hy- 
brid model calculation. Since the abundance of strange 
particles is relatively small they survive the interactions 
in the UrQMD evolution that follows the hydrodynamic 
freeze-out almost without re-thermalization. 

Fig. 1 141 shows the midrapidity yields of protons, pions, 
A's and kaons as a function of the beam energy for central 
Au+Au/Pb+Pb collisions at E lah = 2- 160A GeV. For 
pions, kaons and A's the same trend as for the 47t mul- 
tiplicities is observed. There are less pions produced in 
the hybrid model calculation due to entropy conservation, 
but more strange particles because of the production ac- 
cording to the local thermal equilibrium distributions. 
The proton yield at midrapidity is very similar in both 
calculations while there are less antiprotons produced in 
the hybrid calculation. Also for the strange antiparti- 
cles a reduction of the midrapidity yield in the hybrid 
calculation at higher SPS energies can be seen. 

To explore the kinetics of the system in more detail, 
Fig. [15] shows the rapidity distribution for 7r~ at three 
different energies (E iah = 11,40 and 160 A GeV). The 
general shape of the distribution is very similar in both 
approaches and in line with the experimental data. At 
higher energies even the absolute yields become very close 
to each other in both approaches. 

Fig. [TH] shows the K + rapidity distributions. In this 
case, the yield is higher in the hybrid calculation as al- 
ready discussed above and also the shape of the distri- 
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FIG. 13: (Color online) Excitation function of particle mul- 
tiplicities (47r) in Au+Au/Pb+Pb collisions from _Ei a b = 
2A GeV to = 200 GeV. UrQMD+Hydro (HG) cal- 

culations are depicted with full lines, while UrQMD-2.3 cal- 
culations are depicted with dotted lines. The corresponding 
data from different experiments @, H2, [H [H, [H, 13, HE HE 
[Tool . Hoil . [Tol] are depicted with symbols. 



bution fits very nicely to the experimental data at SPS 
energies. Overall the rapidity distributions seem not to 
be too sensitive to the details of the dynamics for the hot 
and dense stage, but strangeness yields are influenced by 
the local equilibrium assumption. It seems that the local 
equilibrium assumption provides similar strangeness en- 
hancement as previous calculations including additional 
strong color fields [67|, [96| . It is remarkable how well the 
hybrid calculation matches the rapidity spectra at lower 
energies (£7iab = 11A GeV), even though the transport 
calculation provides a slightly better description to the 
experimental data at this energy. One might still con- 
clude from the rapidity spectra that the local equilibrium 
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FIG. 14: (Color online) Excitation function of particle yields 
at midrapidity (\y\ < 0.5) in Au+Au/Pb+Pb collisions from 
Slab = 2A GeV to y/^]7 = 200 GeV. UrQMD+Hydro (HG) 
calculations are depicted with full lines, while UrQMD-2.3 
calculations are depicted with dotted lines. T he c o rrespond- 
ing, data from different experiments HE HE EH GSE GM 
[i05l . QUI, [l07|, [Tol, [Tol, GUI are depicted with symbols. 



is not a good assumption for AGS energies. 



B. Transverse Dynamics 

After the longitudinal dynamics which reflects more 
the stopping power in the initial state we turn now to the 
transverse dynamics of the system. Transverse spectra 
are a promising candidate to be sensitive to the change 
in the underlying dynamics because they emerge from 
the transverse expansion which is mostly dominated by 
the evolution in the hot and dense stage of the reaction. 
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FIG. 15: (Color online) Rapidity spectra of it for cen- 
tral (b < 3.4 fm) Au+Au/Pb+Pb collisions for £i ab = 
11,40 and 160 A GeV. UrQMD+Hydro (HG) calculations are 
depicted with full lines, while UrQMD-2.3 calculations are 
depicted with dotted li nes. The corresponding data from dif- 
ferent experiments [95l . Illl| are depicted with symbols. 



FIG. 16: (Color online) Rapidity spectra of K for cen- 
tral (b < 3.4 fm) Au+Au/Pb+Pb collisions for E^ h = 
11,40 and 160A GeV. UrQMD+Hydro (HG) calculations are 
depicted with full lines, while UrQMD-2.3 calculations are 
depicted with dotted li nes. The corresponding data from dif- 
ferent experiments [95l . Illlfl are depicted with symbols. 



Figs. IT7T [TBI and [PHI display the transverse mass spectra 
for pions, protons and kaons at midrapidity for central 
Au+Au/Pb+Pb reactions at three different beam ener- 
gies. At E\ ah = 11 A GeV (Fig. [TTJ) the differential trans- 
verse mass spectra are very similar for both calculations 
and are in line with the experimental data. 

At = 40,4 GeV (Fig. QJJ) first differences become 
visible. Most notably is the strong flow of protons in 
the hybrid approach, that results in an overestimate of 
protons at high transverse momenta. 

At the highest SPS energy (E iah = 160A GeV, Fig. Q3J) 
all the transverse mass spectra are flatter in the hybrid 
approach. The initial pressure gradients are higher in 
the hydrodynamic calculation due to the hadronic equa- 
tion of state without phase transition. Therefore, the 
matter expands faster in the transverse plane and higher 
transverse masses are reached. At this energy either the 
introduction of a mixed phase (first order phase transi- 
tion) or non-equilibrium effects are necessary to explain 
the experimental data. 

Fig. [201 shows the mean transverse mass excitation 
function for pions. It confirms the observations from the 
differential spectra. Up to 10 GeV beam energy the hy- 
brid model calculation leads to similar results as the de- 
fault UrQMD calculation and is in line with the experi- 
mental data. The mean value of the transverse mass of 
pions is proportional to the temperature of the system 



and very different in the two calculations at higher en- 
ergies. The UrQMD approach shows a softening of the 
equation of state in the region where the phase transition 
is expected because of non-equilibrium effects, while the 
hadron gas hydrodynamic calculation continuously rises 
as a funct ion o f the energy. This behaviour is well known, 
see, e.g., [Till- 
Finally Fig. [2T]shows the mean transverse momenta as 
a function of particle mass for 7r, K, p, A, 3, and f2 parti- 
cles. Here we observe the behaviour known from previous 
hybrid studies, that with increased strangeness, baryons 
accumulate less flow than in a complete hydrodynamic 
approach. This effect can be traced back to the small 
cross sections of multi strange baryons in the hadronic 
cascade, thus showing, that the freeze-out/decoupling 
process proceeds gradually. 



VI. SUMMARY 

We have presented the first fully integrated Boltz- 
mann+hydrodynamics approach to relativistic heavy ion 
reactions. This hybrid approach is based on the Ultra- 
relativistic Quantum Molecular Dynamics (UrQMD) 
transport approach with an intermediate hydrodynam- 
ical evolution for the hot and dense stage of the colli- 
sion. The specific coupling procedure including the ini- 
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FIG. 17: (Color online) Transverse mass spectra of 7r , K + 
and protons at midrapidity (\y\ < 0.5) for central (b < 3.4 fm) 
Au+Au collisions at £ lab = 11 A GeV. UrQMD+Hydro (HG) 
calculations are depicted with full lines, while UrQMD-2.3 
calculations are depicted with dotted l ines. The correspond- 
ing data from the E866 experiment [92l 1 1061 . 1 1 1 it ] are depicted 
with symbols. 



FIG. 19: (Color online) Transverse mass spectra of n ,K + 
and protons at midrapidity (\y\ < 0.5) for central (6 < 3.4 
fm) Pb+Pb collisions for £ lab = 1604 GeV. UrQMD+Hydro 
(HG) calculations are depicted with full lines, while UrQMD- 
2.3 calculations are depicted with dotted lines. Th e corre- 
sponding data from the NA49 experiment [95l . Ill2j ] are de- 
picted with symbols. 
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FIG. 18: (Color online) Transverse mass spectra of 7r , K + 
and protons at midrapidity (\y\ < 0.5) for central (b < 3.4 fm) 
Pb+Pb collisions at E lllb = 40A GeV. UrQMD+Hydro (HG) 
calculations are depicted with full lines, while UrQMD-2.3 
calculations are depicted with dotted lines . Th e correspond- 
ing data from the NA49 experiment [95l . Ill2| are depicted 
with symbols. 



tial conditions and the freeze-out prescription have been 
explained. The event-by-event character of the hybrid 
approach has been emphasized. The present implemen- 
tation allows to compare pure microscopic transport cal- 
culations with hydrodynamic calculations using exactly 
the same initial conditions and frcczc-out procedure. 

The parameter dependences of the model have been 
investigated and the time evolution of different quantities 
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FIG. 20: (Color online) Mean transverse mass excitation func- 
tion of pions at midrapidity (\y\ < 0.5) for central (b < 3.4 
fm) Au+Au/Pb+Pb collisions from £i ab = 2 - 160A GeV. 
UrQMD+Hydro (HG) calculations are depicted with full 
lines, while UrQMD-2.3 calculations are depicted with dot- 
ted lin es. T he corresponding data from different experiments 
0, HE, i.l'l.'i are depicted with symbols. 



were explored. These tests led to the conclusion that the 
choice of the starting time and the freeze-out criterium 
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FIG. 21: (Color online) Mean transverse momentum at 
midrapidity (\y\ < 0.5) as a function of the particle mass 
for pions, kaons, protons, A, 3 and Q in central (b < 3.4 fm) 
Pb+Pb collisions at E lah = 40A GeV. UrQMD+Hydro (HG) 
calculations are depicted with blue squares, while UrQMD-2.3 
calculations are depicted with red circles. 



The effects of the change in the underlying dynam- 
ics - ideal fluid dynamics vs. non-equilibrium transport 
theory - have been explored. The final pion and proton 
multiplicities are lower in the hybrid model calculation 
due to the isentropic hydrodynamic expansion while the 
yields for strange particles are enhanced due to the local 
equilibrium in the hydrodynamic evolution. The results 
of the different calculations for the mean transverse mass 
excitation function, rapidity and transverse mass spec- 
tra for different particle species at three different beam 
energies have been discussed in the context of the avail- 
able data. The transverse expansion of the system is 
much faster in the hybrid model calculation, especially 
at higher energies which leads to differences in the ob- 
servables that are sensitive to the transverse dynamics. 
This finding indicates qualitatively that "new" physical 
effects like, e.g., non-equilibrium effects or a phase tran- 
sition have to be taken into account. 

Forthcoming work will be devoted to the study of dif- 
ferent equations of state and the effect of changes in the 
equation of state on observables. Also in progress are 
calculations within the hybrid model at higher beam en- 
ergies (RHIC and LHC), but therefore specific numerical 
subtleties have to be resolved like, e.g., a dynamical grid 
size for the hydrodynamical evolution. 
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